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ABSTRACT 

Two complex geometry problems are solved using the linearized Euler equations. The impedance mismatch 
method 1 is used to impose the solid surfaces without the need to use a body-fitted grid. The problem is solved in 
the frequency domain to avoid long run times. Although the harmonic assumption eliminates all time dependence, a 
pseudo-time term is added to allow conventional iterative methods to be employed. A Jameson 2 type, Runge-Kutta 
scheme is used to advance the solution in pseudo time. The spatial operator is based on a seven-point, sixth-order finite 
difference. Constant coefficient, sixth-derivative artificial dissipation is used throughout the domain. A buffer zone 
technique employing a complex frequency 3 to damp all waves near the boundaries is used to minimize reflections. The 
results show that the method is capable of capturing the salient features of the scattering, but an excessive number of grid 
points are required to resolve the phenomena in the vicinity of the solid bodies because the wavelength of the acoustics 
is relatively short compared with the size of the bodies. Smoothly transitioning into the immersed boundary condition 
alleviates the difficulties, but a fine mesh is still required. 

INTRODUCTION 

Problems involving complex geometries are commonplace in real world applications. Although considerable ad- 
vances have been made in unstructured, multiblock, and overset grid methods, many of the algorithms that can take 
advantage of these capabilities do not have sufficient accuracy for acoustic problems. Even for high- accuracy methods 
with these capabilities, generating the grid is non-trivial, especially in three-dimensions. For acoustic problems involv- 
ing only scattering, Chung and Morris 1 introduced an alternative that allows a Cartesian mesh to be used for arbitrarily 
complex geometries. All that is required to simulate solid surfaces is to change the density in the regions occupied by 
the bodies. Although very simple to apply, there are some numerical and resolution issues that will be illustrated in the 
example problems. 

Another difficulty in most acoustic problems is the excessive run time required to wash out initial transients and 
obtain a periodic steady- state. For linearized problems that are truly periodic, the time dependence can be factored 
out and the equations solved using algorithms used to accelerate solutions to a steady- state. 4 Unfortunately, there are 
difficulties in applying standard acceleration algorithms to the resulting complex equations. Furthermore, for problems 
with broad frequency content, it is likely to be much more efficient to solve the problem in the time domain with a single 
computation. However, there are some classes of problems where only the behavior of a few important frequencies are 
needed. 

This paper proceeds by describing a numerical method for solving the harmonic form of the linearized Euler equa- 
tions and an implementation of the impedance mismatch method for simulating solid bodies. Example problems with 
two and three cylinders and a periodic source are solved to investigate the merits and disadvantages of the approaches. 
Some conclusions are made on the efficiency of the method as well as possible extensions. 

NUMERICAL METHOD 

The linearized Euler equations for a constant mean flow can be written in the form 
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as reference quantities. The ratio of specific heats is denoted by 7. The source term is given by 

5 = exp — ln(^)~g^4~~ sin ( 87rt ) ( 2 ) 

Equations 1 can be placed in the semi-discrete, compact form 

^ = -7(Q)-D(Q)]. (3) 

Here, T represents the discrete form of all the spatial derivatives and the source term S appearing in equations 1 . A 
seven point, sixth-order operator is used for all spatial derivatives. V denotes sixth-derivative, constant coefficient, 
artificial dissipation. If we assume that all variables have a temporal time dependence of the form exp(— iwt), equation 
3 becomes 

— icjQ exp(— iut) = — [jT(Q) — £>(Q)] exp(— iut). (4) 

Equation 4 can be solved directly, or an iterative technique can be employed. Direct inversion is the most computa- 
tionally efficient approach, but the memory requirement can be excessive. Research is ongoing into parallel versions of 
sparse matrix solvers that may make direct inversion feasible for large problems. However, In the current implementa- 
tion, the equations are solved iteratively by adding a pseudo-time yielding 

^ = - 7(Q) - V(Q) - iwQ] = -K(Q). (5) 

Equation 5 is advanced using a Jameson 2 type, five-stage, Runge-Kutta scheme of the form 

Q (0) = Q ra , 

Q(s) = _ 7(Q(*-D) - P(Q(°>) - , 

1 + a s Lu/\t 

Q n+ 1 = (6) 

where n is the iteration counter and s = 1...5 indicates the stage. The values of a s are [1/4, 1/6, 3/8, 1/2, 1], It 
should be noted that the stability characteristics of the scheme are considerably different in the complex plane, and the 
coefficients used in the Runge-Kutta scheme have been modified appropriately. To accelerate convergence, local time- 
stepping and implicit residual smoothing 5 are employed. Multigrid acceleration is generally less effective with complex 
equations, and is not appropriate for immersed boundary conditions because the effective location of the boundary varies 
with the mesh level. 

Boundary Conditions 

The impedance mismatch method 1 is used to impose solid bodies within the flow. All that is required is to set 
the mean density and pressure to 1/30 of their normal values when within solid bodies. Although one would normally 
expect that using a higher density within solid bodies would yield the desired solution, this problem is difficult to solve 
numerically because spatial derivatives of the acoustic variables must balance the discontinuities in the mean quantities. 
Instead, one can cast the governing equations in the form of 1 where the mean quantities are within the spatial derivatives 
producing smooth fluxes across the solid interfaces. The correct reflected wave is produced, but the solutions within 
solid bodies are fictitious. 

The specified problems are symmetric about y = 0, so symmetry is applied using a reflection condition. All other 
boundaries employ a buffer domain technique 3 to minimize reflections. The frequency in equation 6 is made complex in 
the buffer so that waves are damped. The buffer is four wavelengths thick, and the complex component of the frequency 
varies as 07 $= 2u r (1 — exp(2d 2 ))/ (1 — exp (2)) where d varies from unity at the start of the buffer and is zero at the 
outer boundary. 

RESULTS 

Scattering from Two Cylinders 

Instantaneous pressure contours for the problem with two cylinders is presented in figure 1. In all cases, a uniform, 
Cartesian grid is used for the calculations. The grid extended 25 cylinder diameters in the x direction and 8 diameters in 
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the y direction. The reference cylinder at (-4,0) has a diameter of unity, whereas the cylinder at (4,0) has a diameter of 
0.5. Because the impedance mismatch condition is used to simulate the solid cylinders, there is still a solution within the 
cylinders. However, the values are meaningless. The figure illustrates the complex radiation pattern that results because 
of the interaction of the scattered fields from the two cylinders and the direct radiation from the source. The lack of 
contours near the maximum y boundary in figure 1 reveals that the buffer zone boundary condition is damping out the 
waves. The damped region is not observed in the x direction because the plot has been truncated in that direction. 



Figure 1 : Instantaneous pressure contours for two cylinder problem. 

Figure 2 shows averaged p 2 results along the centerline y — 0. From figure 2(a), it it can be seen that there are some 
differences between the 15 and 30 points-per- wavelength (PPW) solutions. The PPW are measure along a grid line, so 
there is less resolution along the diagonal. Figure 2(b) concentrates on the region around the cylinder at (-4,0). The 
symbols represent the analytic solution. Clearly, neither computational solution is in good agreement with the analytic 
solution in the vicinity of the cylinder. The amplitude of the oscillations on the side of the cylinder closest to the source 
(x > —3.5) are too low, whereas they appear too high in the shadow region (x < —4.5). 



(a) Global view 


(b) Close-up 


Figure 2: Averaged p 2 results along y = 0 for two cylinder problem. 


The numerical algorithm can resolve acoustic wave propagation with considerably less than 15 PPW, so the errors 
are clearly associated with the representation of the cylinders. The actual location of the solid body is somewhat 
ambiguous with the impedance mismatch method because it occurs between grid lines. Therefore, the boundary changes 
as the grid is refined. In the current problem, the acoustic wavelength is considerably smaller than the diameter of the 
cylinders which exacerbates the problem. In addition, although the fluxes are continuous across the interface, the 
solution variables themselves are not. Hence, there is a difficulty in applying artificial dissipation in this region. This 
produces numerical oscillations which can contaminate the solution. Furthermore, the solution in the vicinity of the 
solid bodies is often less accurate than in the far-held. All of these issues contribute to the errors seen in figure 2. 


NAS A/CP— 2004-2 1 2954 


293 


One way to alleviate the numerical errors is to smooth the transition to the lower density within solid bodies. Figure 
3(a) compares the density variation in the vicinity of the cylinder using the standard method which abruptly transitions 
across one grid cell (like a Heaviside function) with that of a 1/2(1 + cos) function that transitions over seven points. 
The results near the cylinder in figure 3(b) shows that the solution with the smooth transition is in considerably better 
agreement with the analytic solution than the one employing a Heaviside function transition. However, the method of 
smoothly varying the density only works well when the distance covered by the transition region is small compared to 
the acoustic wavelength. Furthermore, it can be difficult to determine in which direction to apply the smooth transition 
for a general body. Nonetheless, it can produce considerably better results as seen in figure 4 which presents the pressure 
results on the cylinder surfaces. Because no grid points are actually on the cylinder surfaces, values were interpolated to 
the surface locations. Considering that the exact location of the surface is ambiguous, the results are reasonably good. 




(a) Transition Function 


(b) Solution near cylinder 


Figure 3: Transition function effect on averaged p 2 results along y = 0 for two cylinder problem. The grid has 30 PPW. 




(a) Cylinder at (—4, 0) (b) Cylinder at (4, 0) 

Figure 4: Averaged p 2 results on cylinder surfaces for two cylinder problem. 

Scattering from Three Cylinders 

Instantaneous pressure contours for the problem with three cylinders are presented in figure 5. The grid extended 23 
cylinder diameters in the x direction and 12 diameters in the y direction. The reference cylinder at (-3,0) has a diameter 
of unity, whereas the cylinders at (3,4) and (3,-4) have diameters of 0.75. Considerably more interference is observable 
compared with the two cylinder solution shown in figure 1 . 

Figure 6 shows averaged p 2 results along the centerline y = 0. Although there are some differences between the 
solutions with 15 and 30 PPW, they are not as significant as for the two cylinder problem. The close-up view of a 
solution with a smooth density transition in figure 6(b) does show better agreement with the analytic solution, but it is 
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Figure 5: Instantaneous pressure contours for three cylinder problem. 


only marginally better than the standard solution with a Heaviside type variation. The orientation and size of the smaller 
cylinders relative to the larger one may make the effect of secondary reflections less of an issue along the line y = 0. 




(b) Close-up, 30 PPW 


Figure 6: Averaged p 2 results along y — 0 for three cylinder problem. 


Surface pressure comparisons are presented in figure 7. The solution on the larger cylinder (figure 7(a)) for the case 
with the smooth density transition does agree better with the analytic solution, but the solution on the smaller cylinder 
(figure 7(b)) is considerably more wavy than the solution using the standard method. It is unclear why the solution with 
the smooth density transition exhibits such waviness. Overall, the solutions for the three-cylinder case seem reasonable, 
although the grids are still excessively fine in order to resolve the regions around the cylinders. 

RESOURCES 

The CFL number for the calculations was 2.5. Three million grid points and 6241 iterations were used for the two 
cylinder problem, and 3.8 million points and 6998 iterations for the three-cylinder problem. The residual was forced to 
drop to le-10 in all cases. However, the solution doesn’t vary significantly once the residual drops by several orders of 
magnitude. Hence, one could decrease the number of iterations by a factor of two and obtain the same basic solution. 

A Beowulf cluster of 2.53 GHz Intel Pentium IV computers was used to run the calculations. The two cylinder 
problem took 16 hours on 36 processors, and the three cylinder problem took 25 hours on the same number of processors. 
The code used the Message Passing Interface (MPI) to perform communication over a standard Ethernet interconnect. 
The parallel efficiency is known to be quite low because of a communications bottleneck. 
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(a) Cylinder at (—3, 0) (b) Cylinder at (3, 4) 

Figure 7: Averaged p 2 results on cylinder surfaces for three cylinder problem. The grid has 30 PPW. 

CONCLUSIONS 

The performance of the method used in this work is rather disappointing, although the time to generate the grid 
and make coding changes were almost negligible. The need to have relatively high resolution in the vicinity of solid 
boundaries necessitated a very fine mesh everywhere because of the aim of using a simple Cartesian grid. The relatively 
short wavelength of the acoustics relative to the scattering bodies exacerbated the problem. Furthermore, solving the 
equations in the frequency domain only resulted in a minor improvement in efficiency. Since the grid was uniform, local 
time stepping was irrelevant. Implicit residual smoothing did allow a slightly larger time step to be used, but the im- 
provement was probably less than a factor of 3. A different iterative technique employing some sort of implicitness may 
perform better. However, efficient direct solvers with minimal memory requirements hold the most promise. Further- 
more, the impedance mismatch method does not appear to be the best choice for problems where detailed information 
near the scattering bodies is needed. This is especially true for high-frequency phenomena. If only far-held data for 
lower frequencies is needed, it has already been demonstrated to be useful. However, a finite-element formulation or a 
multiblock method may prove more useful in general. 
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